hilo <- hbb_wku_h_xts
hilo <- data.frame(date=index(hilo), coredata(hilo))
hilo <- hilo[529:54816,]
hilo
# which(hilo$date=="2010-10-23 00:00:00") = index 529
# which(hilo$date=="2016-12-31 23:00:00") = index 54816
length(hilo[,1]) # we are left with 54288 lines of data
[1] 54288
54816 - length(hilo[,1]) # we lost 528 values
[1] 528
# removing columns that we are not using
hilo$date.2 <- NULL # another date column
hilo$date.1 <- NULL # another date column
hilo$BGARFU <- NULL # ?
hilo$cfs <- NULL
hilo$DOmgL <- NULL # dissolved oxygen
#hilo$Doper <- NULL # dissolved oxygen
hilo$PAR1 <- NULL # ?
hilo$pH <- NULL # pH
hilo$NTU <- NULL # a different measurement for turbitity
hilo$DOper10 <- NULL # dissolved oxygen
# colnames(hilo) <- c("Date", "cfs", "RiverFlow-cumec", "LogRiverFlow-cumec", "Chlorophyll-RFU", "Salinity-PPT", "Temp-C", "chlorophyll-calibrator", "Turbidity-NTU")
# does not work ???
hilo
====================================================
length(hilo$logcms[which(is.na(hilo$logcms)==TRUE)]) # 12 NAs
[1] 12
which(is.na(hilo$logcms)==TRUE)
[1] 50509 50510 50511 50512 50513 50514 50515 50516 50517 50518 50519 50520
RiverFlow <- ggplot(hilo, aes(x = date, y = as.numeric(cms))) +
geom_line()
print(RiverFlow + ggtitle("River Flow")+labs(x="Time", y = "River Flow - cubic meters per second"))
length(hilo$ChlRFU[which(is.na(hilo$ChlRFU)==TRUE)]) # 20464 NAs
[1] 20464
which(as.numeric(hilo$ChlRFU)==max(as.numeric(na.omit(hilo$ChlRFU)))) # 15.3 max
[1] 38974
# CHL tells us where in the data set this happened
hilo[38974,]
CHL <- ggplot(hilo, aes(x = date, y = as.numeric(ChlRFU))) +
geom_line()
print(CHL + ggtitle("Chlorophyll ")+labs(x="Time", y = "Chlorophyll - relative fluorescence units (RFU)"))
length(hilo$Corr.NTU[which(is.na(hilo$Corr.NTU)==TRUE)]) #15012 NAs
[1] 15012
which(as.numeric(hilo$Corr.NTU)==max(as.numeric(na.omit(hilo$Corr.NTU)))) # 88.4
[1] 33243
# tells us where in the data set this happened
hilo[33243,]
TURB <- ggplot(hilo,aes(x = date, y = as.numeric(Corr.NTU))) +
geom_line()
print(TURB + ggtitle("Turbidity ")+labs(x="Time", y = "Turbidity - Nephelometric Turbidity Units (NTU)"))
length(hilo$saltppt[which(is.na(hilo$saltppt)==TRUE)]) #11330 NAs
[1] 11330
SALT <- ggplot(hilo, aes(x = date, y = as.numeric(saltppt))) +
geom_line()
print(SALT + ggtitle("Salinity")+labs(x="Time", y = "Salinity - unit parts per thousand (PPT)"))
======================================================== # Histograms FULL DATA SET ## River Flow Histogram
hist(as.numeric(hilo$logcms), main = "Histogram of Log River Flow", xlab = "Log River Flow")
# this looks okay
# VERY skewed
hist(as.numeric(hilo$ChlRFU), main = "Histogram of Chlorophyll", xlab = "Chlorophyll - relative fluorescence units (RFU)")
# this looks better
hist(log(as.numeric(hilo$ChlRFU)), main = "Histogram of Log Chlorophyll", xlab = "Chlorophyll")
NaNs produced
# not sure what happens to units when taking the log
# VERY skewed
hist(as.numeric(hilo$Corr.NTU), main = "Histogram of Turbidity", xlab = "Turbidity - Nephelometric Turbidity Units (NTU)")
# this looks better
hist(log(as.numeric(hilo$Corr.NTU)), main = "Histogram of Log Turbidity", xlab = "Turbidity")
# not sure what happens to units when taking the log
# skewed
hist(as.numeric(hilo$saltppt), main = "Histogram of Salinity", xlab = "Salinity - unit parts per thousand (PPT)")
# this is worst!
hist(log(as.numeric(hilo$saltppt)), main = "Histogram of Log Salinity", xlab = "Salinity")
# not sure what happens to units when taking the log
========================================================= # NEW DATA SET-Modified Data 2013-2015
# start date: 2013-01-01 00:00:00
# end date: 2015-12-31 23:00:00
hilomodified <- hilo[19225:45504,]
length(hilo[,1])-length(hilomodified[,1])
[1] 28008
# lost 28008 entries of data
length(hilomodified[,1])-528 # we are left with 25752 entries of data
[1] 25752
head(hilomodified)
tail(hilomodified)
NA
library(mosaic)
favstats(hilomodified$cms)
Auto-converting character to numeric.
favstats(hilomodified$ChlRFU)
Auto-converting character to numeric.
favstats(hilomodified$Corr.NTU)
Auto-converting character to numeric.
favstats(hilomodified$saltppt)
Auto-converting character to numeric.
favstats(hilomodified$TempC)
Auto-converting character to numeric.
favstats(hilomodified$Doper)
Auto-converting character to numeric.
===================================================== # MODIFIED DATA SET 2013-2015 # Descriptives: Plots
length(hilomodified$logcms[which(is.na(hilomodified$logcms)==TRUE)]) # No NAs, YAY!
[1] 0
which(is.na(hilomodified$logcms)==TRUE)
integer(0)
RiverFlowMod <- ggplot(hilomodified, aes(x = date, y = as.numeric(cms))) +
geom_line()
print(RiverFlowMod + ggtitle("River Flow")+labs(x="Time", y = "River Flow - cubic meters per second"))
sum(is.na(hilomodified$ChlRFU)==TRUE) # 1884 NAs
[1] 1884
#which(is.na(hilomodified$ChlRFU)==TRUE)
# max(as.numeric(na.omit(hilo$ChlRFU))) # 15.3
which(as.numeric(hilomodified$ChlRFU)==15.3)
[1] 19750
hilo[38974,]
CHLMod <- ggplot(hilomodified, aes(x = date, y = as.numeric(ChlRFU))) +
geom_line()
print(CHLMod + ggtitle("Chlorophyll ")+labs(x="Time", y = "Chlorophyll - relative fluorescence units (RFU)"))
length(hilomodified$Corr.NTU[which(is.na(hilomodified$Corr.NTU)==TRUE)]) # 2704 NAs
[1] 2704
# max(as.numeric(na.omit(hilo$Corr.NTU))) # 88.4
which(as.numeric(hilo$Corr.NTU)==88.4)
[1] 33243
hilo[33243,]
TURBMod <- ggplot(hilomodified, aes(x = date, y = as.numeric(Corr.NTU))) +
geom_line()
print(TURBMod + ggtitle("Turbidity ")+labs(x="Time", y = "Turbidity - Nephelometric Turbidity Units (NTU)"))
length(hilomodified$saltppt[which(is.na(hilomodified$saltppt)==TRUE)]) # 2267 NAs
[1] 2267
SALTMod <- ggplot(hilomodified, aes(x = date, y = as.numeric(saltppt))) +
geom_line()
print(SALTMod + ggtitle("Salinity")+labs(x="Time", y = "Salinity - unit parts per thousand (PPT)"))
length(hilomodified$TempC[which(is.na(hilomodified$TempC)==TRUE)]) # 2267 NAs
[1] 1702
TempMod <- ggplot(hilomodified, aes(x = date, y = as.numeric(TempC))) +
geom_line()
print(TempMod + ggtitle("Temperature")+labs(x="Time", y = "Temperature - Celsius"))
length(hilomodified$Doper[which(is.na(hilomodified$Doper)==TRUE)]) # 2267 NAs
TempMod <- ggplot(hilomodified, aes(x = date, y = as.numeric(Doper))) +
geom_line()
print(TempMod + ggtitle("Dissolved Oxygen")+labs(x="Time", y = "Dissolved Oxygen in percent of saturation"))
=================================================== # Histograms MODIFIED
hist(as.numeric(hilomodified$cms), main = "Histogram of Log River Flow", xlab = "Log River Flow", breaks =90, xlim = c(0,100))
# this looks okay
# VERY skewed
hist(as.numeric(hilomodified$ChlRFU), main = "Histogram of Chlorophyll", xlab = "Chlorophyll - relative fluorescence units (RFU)")
# this looks better
hist(log(as.numeric(hilomodified$ChlRFU)), main = "Histogram of Log Chlorophyll", xlab = "Chlorophyll")
NaNs produced
# not sure what happens to units when taking the log
# VERY skewed
hist(as.numeric(hilomodified$Corr.NTU), main = "Histogram of Turbidity", xlab = "Turbidity - Nephelometric Turbidity Units (NTU)")
# this looks better
hist(log(as.numeric(hilomodified$Corr.NTU)), main = "Histogram of Log Turbidity", xlab = "Turbidity")
# not sure what happens to units when taking the log
# skewed
hist(as.numeric(hilomodified$saltppt), main = "Histogram of Salinity", xlab = "Salinity - unit parts per thousand (PPT)")
# this is worst!
hist(log(as.numeric(hilomodified$saltppt)), main = "Histogram of Log Salinity", xlab = "Salinity")
# not sure what happens to units when taking the log
=============================================================
It is hard to see what is going on
AllYears <- ggplot(hilomodified, aes(x = date, y = as.numeric(logcms))) +
geom_line()+
geom_line(aes(y = as.numeric(saltppt)), color = "darkred") +
geom_line(aes(y = as.numeric(Corr.NTU)), color="darkgreen") +
geom_line(aes(y=as.numeric(ChlRFU)),color="blue")
print(AllYears + ggtitle("Hilo Bay")+labs(x="Time", y = "River Flow - cubic meters per second"))
============================================================== # Descriptives by Storm We are picking one storm from each year. We can indicate a storm has occurred by the extreme events in the river flow data. We will not use the log (which is logbase10) in order to see the extreme events When salinity is below 35 this also indicates a storm has occurred.
We will break the data set by year to find the most extreme event for each year.
hilo2013 <- hilomodified[1:8760,]
hilo2014 <- hilomodified[8761:17520,]
hilo2015 <- hilomodified[17521:length(hilomodified[,1]),]
# max(as.numeric(hilo2013$logcms)) # this is log base 10
# This is NOT the natural log
max(as.numeric(hilo2013$cms)) # 207.186
[1] 207.186
which(as.numeric(hilo2013$cms) == max(as.numeric(hilo2013$cms)))
[1] 3550
hilo2013[3550,]
# all values for this storm are NA
plot2013 <- ggplot(hilo2013, aes(x = date, y = as.numeric(cms))) +
geom_line(color="black")+
geom_line(aes(y = as.numeric(saltppt)), color = "darkred") +
geom_line(aes(y = as.numeric(Corr.NTU)), color="darkgreen") +
geom_line(aes(y=as.numeric(ChlRFU)),color="blue")
print(plot2013 + ggtitle("2013")+labs(x="Time", y = "River Flow - cubic meters per second"))
# R2013.1 <- ggplot(hilo2013[1:(length(hilo2013[,1])/2),], aes(x = date, y = as.numeric(cms))) +
# geom_line(color="black")+
# geom_line(aes(y = as.numeric(saltppt)), color = "darkred") +
# geom_line(aes(y = as.numeric(Corr.NTU)), color="darkgreen") +
# geom_line(aes(y=as.numeric(ChlRFU)),color="blue")
# print(R2013.1 + ggtitle("River Flow")+labs(x="Time", y = "River Flow - cubic meters per second"))
# R2013.1+ylim(0,40)
#
#
# R2013.2 <- ggplot(hilo2013[(length(hilo2013[,1])/2):length(hilo2013[,1]),], aes(x = date, y = as.numeric(cms))) +
# geom_line(color="black")+
# geom_line(aes(y = as.numeric(saltppt)), color = "darkred") +
# geom_line(aes(y = as.numeric(Corr.NTU)), color="darkgreen") +
# geom_line(aes(y=as.numeric(ChlRFU)),color="blue")
# print(R2013.2 + ggtitle("River Flow")+labs(x="Time", y = "River Flow - cubic meters per second"))
# R2013.2+ylim(0,40)
# Max river flow in the overall data set
max(as.numeric(hilo2014$cms))
[1] 472.4939
max(as.numeric(hilomodified$cms))
[1] 472.4939
which(as.numeric(hilo2014$cms) == max(as.numeric(hilo2014$cms)))
[1] 5263
hilo2014[5263,]
R2014 <- ggplot(hilo2014, aes(x = date, y = as.numeric(logcms))) +
geom_line(color="black")+
geom_line(aes(y = as.numeric(saltppt)), color = "darkred") +
geom_line(aes(y = as.numeric(Corr.NTU)), color="darkgreen") +
geom_line(aes(y=as.numeric(ChlRFU)),color="blue")
print(R2014 + ggtitle("2014")+labs(x="Time", y = "River Flow - cubic meters per second"))
max(as.numeric(hilo2015$cms))
[1] 232.2449
which(as.numeric(hilo2015$cms) == max(as.numeric(hilo2015$cms)))
[1] 6475
hilo2015[6475,]
R2015 <- ggplot(hilo2015, aes(x = date, y = as.numeric(logcms))) +
geom_line(color="black")+
geom_line(aes(y = as.numeric(saltppt)), color = "darkred") +
geom_line(aes(y = as.numeric(Corr.NTU)), color="darkgreen") +
geom_line(aes(y=as.numeric(ChlRFU)),color="blue")
print(R2015 + ggtitle("2015")+labs(x="Time", y = "River Flow - cubic meters per second"))
outcomes <- as.numeric(hilomodified$saltppt)<25
index.lt25 <- which(outcomes==TRUE)
poss.storms <- hilomodified[index.lt25,]
poss.storms
outcomes.35 <- as.numeric(hilomodified$saltppt)<35
index.lt35 <- which(outcomes==TRUE)
index.lt25
[1] 159 160 161 162 163 164 209 210 213 218 225 227 228
[14] 230 231 232 234 759 765 766 815 816 894 895 899 900
[27] 901 902 903 904 905 906 907 908 1112 1113 1114 1115 1116
[40] 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1131 1132 1133
[53] 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146
[66] 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159
[79] 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172
[92] 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185
[105] 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198
[118] 1200 1201 1202 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213
[131] 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226
[144] 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239
[157] 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252
[170] 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265
[183] 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278
[196] 1279 1283 1284 1285 1286 1288 1289 1293 1294 1295 1296 1297 1298
[209] 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311
[222] 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324
[235] 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337
[248] 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350
[261] 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363
[274] 1364 1365 1366 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377
[287] 1378 1379 1380 1381 1382 1383 1384 1391 1392 1393 1410 1411 1413
[300] 1414 1415 1416 1417 1418 1419 1420 1421 1424 1425 1427 1428 1431
[313] 1432 1433 1434 1435 1436 1437 1438 1439 1440 1442 1443 1444 1445
[326] 1446 1447 1450 1451 1452 1453 1454 1458 1478 1479 1480 1481 1482
[339] 1483 1484 1490 1515 2137 3485 3487 3488 3489 3490 3491 3492 3495
[352] 3496 3497 3498 3499 3500 3501 3502 3503 3508 3509 3513 3518 3519
[365] 5312 5313 5317 5318 5319 5320 5321 5322 5323 5325 5326 5327 5328
[378] 5329 5330 5331 5332 5334 5335 5336 5337 5338 5339 5340 5341 5342
[391] 5343 5344 5345 5346 5347 5348 5349 5350 5368 5474 5475 5485 5486
[404] 5487 5488 5489 5491 5492 5509 5756 7545 8672 8673 8674 8675 8676
[417] 8677 8678 8679 8680 8681 8682 8683 8713 8714 8715 8716 8717 8718
[430] 8719 8720 8721 8722 8723 8724 8725 8726 8727 8728 8729 8730 8731
[443] 8732 8733 8734 8735 8736 8737 8738 8739 8740 8741 8742 8743 8744
[456] 8745 8746 8747 8748 8749 8750 8751 8752 8753 8754 8755 8756 8757
[469] 8758 8759 8760 8761 8762 8763 8764 8765 8766 8767 8768 8769 8770
[482] 8771 8772 8773 8774 8775 8776 8835 8836 8837 8839 8840 8841 8842
[495] 8843 8844 8845 8846 8847 8848 8849 8855 8856 8857 8858 8859 8862
[508] 8863 8864 8865 8876 8927 8932 8933 8945 10394 10395 10396 10397 10398
[521] 10399 10400 10401 10402 10403 10404 10405 10406 10407 10408 10410 10411 10413
[534] 10414 10421 10422 10425 10427 10428 10429 10430 10431 10432 10433 10434 10435
[547] 10436 10437 10438 10441 10442 10443 10444 10445 10446 10450 10538 10539 10540
[560] 10541 10544 10545 10546 10547 10552 10553 10554 10555 10556 10557 10558 10559
[573] 10570 10644 10647 10652 10656 10657 10659 10664 10665 10666 10667 10668 10669
[586] 10670 10671 10672 10673 10674 10675 10676 10677 10678 10679 10680 10681 10682
[599] 10683 10684 10685 10686 10687 10688 10689 10690 10691 10692 10693 10694 10695
[612] 10696 10697 10698 10699 10700 10701 10702 10703 10704 10705 10706 10707 10708
[625] 10709 10710 10711 10712 10713 10714 10715 10716 10717 10718 10719 10720 10721
[638] 10722 10723 10724 10725 10726 10727 10728 10729 10730 10731 10732 10733 10734
[651] 10735 10736 10737 10738 10739 10740 10741 10742 10743 10744 10745 10746 10747
[664] 10748 10749 10750 10751 10752 10753 10754 10755 10756 10757 10758 10759 10760
[677] 10761 10762 10763 10764 10765 10766 10767 10768 10769 10770 10771 10772 10773
[690] 10774 10775 10776 10777 10778 10779 10780 10781 10782 10783 10784 10785 10786
[703] 10787 10788 10789 10790 10791 10792 10793 10794 10795 10796 10797 10798 10799
[716] 10800 10801 10802 10803 10804 10805 10806 10807 10808 10809 10810 10811 10812
[729] 10813 10814 10815 10816 10817 10818 10819 10820 10830 10831 10832 10833 10919
[742] 10920 10921 10922 10923 10924 10925 10926 10927 10928 10929 10930 10931 10932
[755] 10933 10934 10935 10936 10937 10938 10939 10940 10941 10942 10943 10944 10945
[768] 10946 10947 10948 10949 10950 10951 10952 10953 10954 10955 10956 10957 10958
[781] 10959 10960 10964 10965 10966 10967 10968 10969 10970 10971 10972 10973 10974
[794] 10975 10976 10977 10978 10979 10980 10981 10982 10983 10984 10985 10986 10987
[807] 10988 10989 10990 10991 10992 10993 10994 10995 10997 10998 11000 11001 11002
[820] 11003 11004 11005 11006 11007 11008 11009 11010 11011 11012 11019 11020 11021
[833] 11022 11023 11024 11025 11026 11027 11028 11030 11031 11032 11033 11065 11077
[846] 11086 11093 11094 11095 11096 11097 11098 11099 11100 11101 11102 11103 11104
[859] 11105 11106 11107 11108 11109 11110 11111 11112 11113 11114 11115 11116 11117
[872] 11118 11119 11120 11122 11123 11124 11125 11126 11127 11128 11129 11130 11131
[885] 11132 11133 11134 11135 11136 11139 11140 11149 11150 11151 11204 11205 11210
[898] 11211 11212 11213 11214 11215 11216 11217 11218 11219 11220 11221 11222 11223
[911] 11224 11225 11226 11227 11228 11229 11230 11231 11232 11233 11234 11235 11236
[924] 11237 11238 11239 11240 11241 11242 11243 11244 11245 11246 11247 11248 11249
[937] 11250 11251 11252 11253 11254 11255 11256 11257 11258 11259 11260 11261 11262
[950] 11263 11264 11265 11266 11267 11268 11269 11270 11271 11272 11273 11274 11275
[963] 11276 11277 11278 11279 11280 11281 11282 11283 11284 11285 11286 11287 11288
[976] 11289 11290 11291 11292 11293 11294 11295 11296 11297 11298 11299 11300 11301
[989] 11302 11303 11304 11305 11306 11307 11308 11309 11310 11311 11312 11313
[ reached getOption("max.print") -- omitted 2851 entries ]
poss.storms35 <- hilomodified[index.lt25,]
poss.storms35
hilomodified$date[16]
[1] "2013-01-01 15:00:00 UTC"
storm.1.7.13 <- ggplot(hilomodified[145:168,], aes(x = date, y = as.numeric(cms))) +
geom_line(color="black")+
geom_line(aes(y = as.numeric(saltppt)), color = "lightsteelblue1") +
geom_line(aes(y = as.numeric(Corr.NTU)), color="grey69") +
geom_line(aes(y=as.numeric(ChlRFU)),color="khaki")
print(storm.1.7.13 + ggtitle("Storm 1/7/13")+labs(x="Time"))
RiverFlow1 <- ggplot(hilomodified[1:100,], aes(x = date, y = as.numeric(cms))) +
geom_line()
print(RiverFlow + ggtitle("River Flow")+labs(x="Time", y = "River Flow - cubic meters per second"))
=======
getwd()
[1] "/Users/odalysbar/Documents/MA321-SP21/HiloBay-Adolf/HiloBay_MA321_Spring21"
# write.csv(hilomodified,file="HiloBayNEW13to15.csv", row.names = FALSE)
# any value greater than or equal to 10 cms indicates a storm
pos.neg <- as.numeric(hilomodified$cms) - 10
start <- vector()
end <- vector()
for(i in 1:length(pos.neg)){
if(isTRUE(pos.neg[i] < 0 && pos.neg[i+1] > 0)){
# if goes from + to - then this is the beginning of a storm
start[i] <- i
}else if(isTRUE(pos.neg[i] > 0 && pos.neg[i+1] < 0)){
# if goes from - to + then this is the end of the storm
end[i] <- i
}
}
start<-start[!is.na(start)] # removing NA
start
[1] 141 212 891 1097 1106 3480 3499 3545 5309 7537 8713 8729 10388
[14] 10390 10411 10727 10737 10843 10917 11216 11340 11508 11646 12581 13406 13458
[27] 13551 13655 13866 14010 15233 15287 15722 16054 16727 19777 19802 19920 20066
[40] 20386 20622 20662 21143 21163 22249 22656 22762 22911 23008 23015 23060 23116
[53] 23201 23395 23640 23703 23740 23757 23774 23807 23947 24141 24312 24437 24938
[66] 25005 25246 25474 25492 25972 26107 26187
length(start) # make sure both start and end are the same length
[1] 72
end<-end[!is.na(end)]
end
[1] 162 226 903 1102 1398 3495 3500 3572 5338 7547 8718 8765 10389
[14] 10404 10412 10732 10752 10845 11044 11232 11346 11525 11695 12596 13445 13514
[27] 13648 13670 13888 14111 15235 15288 15750 16148 16749 19790 19811 19938 20087
[40] 20394 20650 20719 21153 21169 22257 22675 22771 22972 23013 23020 23070 23130
[53] 23215 23403 23670 23719 23742 23761 23781 23899 24082 24168 24325 24510 24945
[66] 25233 25328 25482 25499 26093 26117 26214
length(end)
[1] 72
library(ggplot2)
for(i in 1:length(start)){
# the point start[i] is the beginning of a storm
# the point end[i] is the end of the same storm
# we are adding 24 hours to each start and end time
p <- ggplot(hilomodified[(start[i]-24):(end[i]+24),],aes(x = date, y = as.numeric(cms))) +
geom_line(color="black")
p <- p + geom_line(aes(y = range01(as.numeric(Corr.NTU))), color="green")
p <- p + geom_line(aes(y= range01(as.numeric(ChlRFU))),color="blue")
p <- p + scale_y_continuous(name = "First Axis", sec.axis = sec_axis(trans = ~.*.05, name="Second Axis")
)
# geom_line(aes(y = range01(as.numeric(Corr.NTU))), color="green") +
# geom_line(aes(y= range01(as.numeric(ChlRFU))),color="blue") +
# geom_line(aes(y = range01(as.numeric(saltppt))),color="red") +
# geom_line(aes(y = range01(as.numeric(TempC))),color="purple") +
# geom_line(aes(y = range01(as.numeric(Doper))),color="orange")+
p <- p + ggtitle("Storms")+labs(x="Time")
print(p)
}
Error in range01(as.numeric(Corr.NTU)) :
could not find function "range01"
range01 <- function(x, na.rm = TRUE) {
# version 1.1 (9 Aug 2013)
(x - min(x, na.rm = na.rm)) / (max(x, na.rm = na.rm) - min(x, na.rm = na.rm))
}
num <- c(1:length(start))
graphics.off()
for(i in 1:length(start)){
#for(i in c(1)){
jpeg(sprintf("storm-%s.jpeg", num[i]),width=2000, height = 1000)
index <- c(start[i]:(end[i]+24))
date <- c(hilomodified$date[index])
enddate <- hilomodified$date[end[i]+24]-10
endcms <- as.numeric(hilomodified$cms[end[i]+24])
endchl <- as.numeric(hilomodified$ChlRFU[end[i]+24])
endturb <- as.numeric(hilomodified$Corr.NTU[end[i]+24])
endsalt <- as.numeric(hilomodified$saltppt[end[i]+24])
endtemp <- as.numeric(hilomodified$TempC[end[i]+24])
endDO <- as.numeric(hilomodified$Doper[end[i]+24])
par(mfrow=c(6,1))
# RIVER FLOW
plot(x = date,y = as.numeric(hilomodified$cms[index]), type = "l",lwd=2, col="black", xlab="Time", ylab="River Flow cms", cex.lab = 1.5) # Create first plot for River Flow
text(x= enddate, y = endcms, col = "black", adj=1, cex=3, label= "River Flow")
# CHLOROPHYLL
if(sum(is.na(as.numeric(hilomodified$ChlRFU[index])))==length(hilomodified$ChlRFU[index])){
plot(x = 0 ,y = 0, type = "l",lwd=2, col="blue", xlab= "Time", ylab = "Chlorophyll - relative fluorescence units (RFU)", cex.lab = 1.5)
text(x = 0, y = 0, col = "blue", adj=1, cex=3, label= "Chlorophyll")
}else{
plot(x = date ,y = as.numeric(hilomodified$ChlRFU[index]), type = "l",lwd=2, col="blue", xlab= "Time", ylab = "Chlorophyll - relative fluorescence units (RFU)", cex.lab = 1.5)
text(x = enddate, y = endchl, col = "blue", adj=1, cex=3, label= "Chlorophyll")
}
# TURBIDITY
if(sum(is.na(as.numeric(hilomodified$Corr.NTU[index])))==length(hilomodified$Corr.NTU[index])){
plot(x = 0 ,y = 0, type = "l",lwd=2, col="red", xlab= "Time", ylab = "Turbidity - Nephelometric Turbidity Units (NTU)", cex.lab = 1.5)
text(x = 0, y = 0, col = "red", adj=1, cex=3, label= "Turbidity")
}else{
plot(x = date ,y = as.numeric(hilomodified$Corr.NTU[index]), type = "l",lwd=2, col="red", xlab= "Time", ylab = "Turbidity - Nephelometric Turbidity Units (NTU)", cex.lab = 1.5)
text(x = enddate, y = endturb, col = "red", adj=1, cex=3, label= "Turbidity")
}
# SALINITY
if(sum(is.na(as.numeric(hilomodified$saltppt[index])))==length(hilomodified$saltppt[index])){
plot(x = 0 ,y = 0, type = "l",lwd=2, col="green", xlab= "Time", ylab = "Salinity - unit parts per thousand (PPT)", cex.lab = 1.5)
text(x = 0, y = 0, col = "green", adj=1, cex=3, label= "Salinity")
}else{
plot(x = date ,y = as.numeric(hilomodified$saltppt[index]), type = "l",lwd=2, col="green", xlab= "Time", ylab = "Salinity - unit parts per thousand (PPT)", cex.lab = 1.5)
text(x = enddate, y = endsalt, col = "green", adj=1, cex=3, label= "Salinity")
}
# TEMPERATURE
if(sum(is.na(as.numeric(hilomodified$TempC[index])))==length(hilomodified$TempC[index])){
plot(x = 0 ,y = 0, type = "l",lwd=2, col="orange", xlab= "Time", ylab = "Temperature - Celsius", cex.lab = 1.5)
text(x = 0, y = 0, col = "orange", adj=1, cex=3, label= "Temperature")
}else{
plot(x = date ,y = as.numeric(hilomodified$TempC[index]), type = "l",lwd=2, col="orange", xlab= "Time", ylab = "Temperature - Celsius", cex.lab = 1.5)
text(x = enddate, y = endtemp, col = "orange", adj=1, cex=3, label= "Temperature")
}
# DISSOLVED OXYGEN
if(sum(is.na(as.numeric(hilomodified$Doper[index])))==length(hilomodified$Doper[index])){
plot(x = 0 ,y = 0, type = "l",lwd=2, col="purple", xlab= "Time", ylab = "Dissolved Oxygen in percent of saturation", cex.lab = 1.5)
text(x = 0, y = 0, col = "purple", adj=1, cex=3, label= "Dissolved Oxygen")
}else{
plot(x = date ,y = as.numeric(hilomodified$Doper[index]), type = "l",lwd=2, col="purple", xlab= "Time", ylab = "Dissolved Oxygen in percent of saturation", cex.lab = 1.5)
text(x = enddate, y = endDO, col = "purple", adj=1, cex=3, label= "Dissolved Oxygen")
}
dev.off()
}
TRUE+TRUE
[1] 2
sprintf("storm-%s",1 )
[1] "storm-1"